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C " 3 ' Abstract. A framework for the analysis of stochastic models of chemical systems for which the 

f^^ ' deterministic mean-field description is undergoing a saddle-node infinite period (SNIPER) bifurcation 

^\1 ^ is presented. Such a bifurcation occurs for example in the modelling of cell-cycle regulation. It is 

, shown that the stochastic system possesses oscillatory solutions even for parameter values for which 

K 1 the mean-field model does not oscillate. The dependence of the mean period of these oscillations on 

3 ' the parameters of the model (kinetic rate constants) and the size of the system (number of molecules 

^"5 ' present) is studied. Our approach is based on the chemical Fokker-Planck equation. To get some 

r/^ insights into advantages and disadvantages of the method, a simple one-dimensional chemical switch 

, is first analyzed, before the chemical SNIPER problem is studied in detail. First, results obtained by 

I 1 1 solving the Fokker-Planck equation numerically are presented. Then an asymptotic analysis of the 

r^ ' Fokker-Planck equation is used to derive explicit formulae for the period of oscillation as a function 

O l' of the rate constants and as a function of the system size. 

^ , Key words, stochastic bifurcations, chemical Fokker-Planck equation 
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C/3 ' 1. Introduction. Bifurcation diagrams are often used for the analysis of deter- 

, y . ministic models of chemical systems. In recent years, they have also been applied 

^ ' to models of biological (biochemical) systems. For example, the cell-cycle model of 

Tyson et al. [13] is a system of ordinary differential equations (ODEs) which describes 
Q,' the time evolution of concentrations of biochemical species involved in cell-cycle reg- 

ulation. The dependence of the qualitative behaviour on the parameters of the model 
is studied using standard bifurcation analysis of ODEs [21 [3] . Biological systems are 
^ , intrinsically noisy because of the low copy numbers of some of the biochemical species 

QiQ ' involved. In such a case, one has to use stochastic simulation algorithms (SSAs) to 

0> , include the intrinsic noise in the modelling [5]. The SSA models use the same set of 

'^ ' parameters as the ODE models. However, even if we choose the same values of the 

^. , rate constants for both the SSA model and the ODE model, the behaviour of the 

C^^ ' system can differ qualitatively [1]. For example, the transition from the G2 phase to 

^^ , mitosis in the model of Tyson et al. is governed by a SNIPER (saddle-node infinite 

period) bifurcatioiu [13] . In the ODE setting, a SNIPER bifurcation occurs whenever 
a saddle and a stable node collapse into a single stationary point on a closed orbit 
(as a bifurcation parameter is varied). In particular, the limit cycle is born with 
infinite period at the bifurcation point. If we use the SSA model instead of ODEs, 
^ ' the bifurcation behaviour will be altered. As we will see in Section [21 intrinsic noise 

causes oscillations with finite average period at the deterministic bifurcation point. 
Moreover, the stochastic system oscillates even for the parameter values for which 
the deterministic model does not have a periodic solution. Clearly, there is a need to 
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understand the changes in the bifurcation behaviour if a modeller decides to use SSAs 
instead of ODEs. In this paper, we focus on the SNIPER bifurcation. However, the 
approach presented can be applied to more general chemical systems. 

In Section [2l we introduce a simple chemical system for which the ODE model 
undergoes a SNIPER bifurcation. We use this illustrative example to motivate the 
phenomena which we want to study. The study of stochastic chemical systems in 
the neighbourhood of determistic bifurcation points will be done using the chemical 
Fokker-Planck equation [7]. The analysis of the Fokker-Planck equation is much 
easier in one dimension. Thus, we start with an analysis of a simple one-dimensional 
chemical switch in Section [S] Using the one-dimensional setting, we can show some 
advantages and disadvantages of our approach without going into technicalities. In 
Section 21 we present a computer assisted analysis of the SSAs by exploiting the 
chemical Fokker-Planck equation, using the illustrative model from Section [2 In 
Section [51 we provide core analytical results. We study the dependence of the period 
of oscillation on the parameters of the model, namely kinetic rate constants and the 
size of the system (number of molecules present). We derive analytical formulae 
for the period of oscillation using the asymptotic expansion of the two-dimensional 
Fokker-Planck equation. Finally, in Section [6l we present our conclusions. 

2. A chemical system undergoing a SNIPER bifurcation. We consider 
two chemical species X and y in a well-mixed reactor of volume V which are subject 
to the following set of seven chemical reactions 



^ y ^ X ^ 0, 2X ^ 3X, (2.1) 



^4d 



X + Y '^ X + 2Y, 2X + Y ^ 2X. (2.2) 

The first reaction is the production of the chemical Y from a source with constant 
rate kid per unit of volume, i.e. the units of kid are sec~^mni~^, see (J2.5I) below. The 
second reaction is the conversion of y to X with rate constant k2d and the third re- 
action is the degradation of X with rate constant k^d- The remaining reactions are of 
second-order or third-order. They were chosen so that the mean-field description cor- 
responding to (|2.ip - (|2.2p undergoes a SNIPER bifurcation as shown below. Clearly, 
there exist many other chemical models with a SNIPER bifurcation, including the 
more realistic model of the cell-cycle regulation |13J. The advantage of the model 
()2.ip -- (|2.2p is that it involves only two chemical species X and Y, making the visu- 
alization of our results clearer: two-dimensional phase planes are much easier to plot 
than the phase spaces of high-dimensional models. The generalization of our results 
to models involving more than two chemical species will be discussed in Section [S] 

Let X = X{t) and Y = Y{t) be the number of molecules of the chemical species 
X and Y, respectively. The concentration of A" (resp. Y) will be denoted by x = X/V 
(resp. y = Y/V). If we have enough molecules of X and Y in the system, we often 
describe the time evolution of x and y by the mean-field ODE model. Using (|2.ip - 
p.2p . this can be written as 

-rr ^ k2dV - k^dx +kidX -k^dX, (2.3) 

at 

-TT ^ -k-jdx y + k^dxy ~ k2dy + kid. (2.4) 

at 
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Fig. 2.1. Nullclines of the ODE system \2.3\ - ^2^ f or kid = 12 (left panel) and kid = 13 
(right panel). The other parameter values are given by 1 12.51 1. The steady states are denoted by blue 
dots. Illustrative trajectories which start close to the steady states are plotted as thin black lines. 



We choose the values of the rate constants as follows: 

kid — 12 [sec~^mni~'^], k2d — 1 [sec~^], k^d — 33 [sec~^], k^d — H [sec" 
ksd = 1 [sec""^ mm®], k^d — 0.6 [sec~^ mm'^], kjd — 0.13 [sec~^ mm®], 



(2.5) 



where we have included the units of each rate constant to emphasize the dependence 
of each rate constant on the volume. From now on, we will treat all rate constants as 
numbers, dropping the corresponding units, to simplify our notation. The nullclines 
of the ODE system (|2.3p - (|2.4p are plotted in Figure [CT aV The x-nullcline (resp. 
2/-nullcline) is given by dx/dt — (resp. dy/dt — 0). The nullclines intersect at three 
steady states which are denoted as SN (stable node). Saddle and UN (unstable node). 
We also plot illustrative trajectories which start close to each steady state (thin black 
lines). As the parameter values change, the stable node (SN) and the saddle can 
approach each other. If they coalesce into one point (with the dominant eigenvalue 
equal to 0), a limit cycle with infinite period appears. Shifting the nullclines further 
apart, we obtain a dynamical system with periodic solutions. This is demonstrated 
in Figure I2.ir b) where we plot the nullclines of (|2.3[) - (|2.4[) for kid — 13, with the 
other parameter values given by (|2.5p . An illustrative trajectory (thin black line) 
converges to the limit cycle. We see that there is a SNIPER (saddle-node infinite 
period) bifurcation as the parameter kid is increased from 12 to 13. 

The main goal of this paper is to understand and analyse changes in the behaviour 
of chemical systems when deterministic ODE models are replaced by SSAs, i.e. when 
the intrinsic noise is taken into account. The stochastic model of the chemical sys- 
tem (|2.ip - (l2.2p is given by the Gillespie SSA |6j which is equivalent to solving the 
corresponding chemical master equation - see Appendix [X] We denote the reaction 
with the rate constant kid, i = 1,2, ... ,7, as the reaction Ri. Note that the reversible 
reaction in (|2.1[) is considered as two separate chemical reactions. To use the Gille- 
spie SSA, we have to specify the propensity function ai{x,y), i — 1,2, ... , 7, of each 
chemical reaction in (j2.ip - ()2.2p . The propensity function is defined so that ai{x, y) di 
is the probability that, given X{t) — x and Y{t) — y, one Ri reaction will occur in the 
next infinitesimal time interval [t, t + At). For example, kid is the rate of production 
of Y molecules per unit of volume. Thus, fci = kidV is the rate of production of Y 
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Fig. 2.2. (a) T/ie time evolution of X given by the stochastic (blue line) and deterministic (red 
line) models of the chemical system i2.1t - i2.2t . The rate constants are given by i2.5t and V = AO. 
(b) The same trajectory plotted in the X-Y plane. 



molecules per the whole reactor of volume V . Consequently, the probability that one 
Y molecule is produced in the infinitesimally small time interval [t, t + dt) is equal 
to kidt, i.e. the propensity function of the reaction Ri is ai{x,y) = fci = ki^V. 
To specify other propensity functions, we first scale the rate constant kid with the 
appropriate power of the volume V, namely we define 



ki = kidV, k2 = k2d, kz = k 



3d, 



kid 

V 



k^d 



k&d 



«7d 



ki = ^^, fcs = -^, h = — , kj^ —. (2.6) 



Then the propensity function of the reaction i?i, i = 1, 2, . . . , 7, is given as the product 
of the scaled rate constant and numbers of available reactant molecules, namely 

ai(a;,y) = /ci, a2{x,y) ^k2V, a3{x,y) ^ ksx, a4{x,y) ^ k4x{x - 1), 



a5{x,y) = k5x{x-l){x-2), ae{x,y) = kexy, a7{x,y) ^ kjx{x - l)y. (2.7) 

Note that the propensity functions of the reactions i?4 and R^ are proportional to the 
number x{x — l)/2 of available pairs of X molecules, and not to x^; similarly, as is 
proportional to the number x{x — l){x ~ 2)/6 of available triplets of X molecules and 
not to x^ . For a general discussion of the propensity functions see [5] . 

Using (|2.7p and the Gillespie SSA, we can simulate the stochastic trajectories 
of (|2.ip - (|2.2p . In Figure I^T^ a). we compare the time evolution of X given by the 
stochastic (Gillespie SSA) and deterministic (the ODE system (I2.3p ~ (l2.4p ) models. 
We use the same parameter values ()2.5p and the same initial condition [X, Y] = 
[40,480] for both, stochastic and deterministic, simulations. The reactor volume is 
V — 40. In Figure [2?2r b). we plot both trajectories in the x-y plane. We see that the 
solution of the deterministic equations converges to a steady state while the stochastic 
model has oscillatory solutions. In Figure 12.31 we present results obtained by the 
stochastic simulation of (|2.1[) - (|2.2[) for kid — 13. The other parameters are the same 
as in Figure [2T2l We see that both the stochastic and deterministic models oscillate 
for kid — 13. 

An important characteristic of oscillating systems is their period of oscillation. 
This is a well-defined number for the deterministic ODE model, but in the stochastic 
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Fig. 2.3. (a) The time evolution of X given by the stochastic (blue line) and deterministic (red 
line) models of the chemical system i2.1t - i2.2t for fei^j = 13. The other parameters are chosen as 
in Figure [2. 2\ (b) The same trajectory plotted in the X-Y plane. 



case the periods of individual oscillations vary. Thus, in the stochastic model we are 
interested in the mean period of oscillation (averaged over many periods). The mean 
period of oscillation is plotted in Figure I2.4r a) as a function of the rate constant kid 
(blue circles). It was computed as an average over 10,000 periods for each of the 
presented values of kid- The starting/finishing point of every period was defined as 
the time when the simulation enters the half-plane X > 200. More precisely, we start 
the computation of each period whenever the trajectory enters the area X > 200. 
Then we wait until X < 60 before we test the condition X > 200 again. The extra 
condition X < 60 guarantees that possible small random fluctuations around the 
point X = 200 are not counted as two or more periods. The period of oscillation of 
the ODE model (|2.3p - (l2.4p is plotted in Figure l^^ a) as the red line. It asymptotes 
to inflnity as {kid — Kd)^^^^ for kid -^ K^ where Kd = 12.2 is the bifurcation value 
of the parameter kid- 
In the limit V ^ oo (which is the so-called thermodynamic limit JJ), the stochas- 
tic description converges to the ODE model (2.3)-(2.4), that is, the probability dis- 
tributions become Dirac-like and their averages converge to the solution of the ODEs 
(2.3)-(2.4) for V -^ oo. The dependence of the period of oscillation on the volume V 
is shown in Figure [^TlT b') where we fix kid equal to the bifurcation value Kd = 12.2 
of the ODE model and we vary the volume V . The other parameter values are given 
by (|2.5p . Since kid — Kd, the period of oscillation of the ODE model is infinity. 
In Figure I2.4r b) , we see that the period of oscillation of the stochastic model is an 
increasing function of V- It approaches the period of oscillation of the ODE model 
(infinity) as y — > cxd. 

The estimates of the period of oscillation (blue circles in Figure 12. 4p were ob- 
tained as averages over 10,000 periods of the Gillespie SSA. Such an approach is 
computationally intensive. The goal of this paper is to show that we can obtain the 
same information by solving and analyzing the chemical Fokker-Planck equation. In 
Section [H we present results obtained by solving the Fokker-Planck equation using 
a suitable finite clement method. In Section [5l we use the asymptotic analysis of 
the Fokker-Planck equation to derive explicit formulae for the period of oscillation, 
as a function of the rate constants and as a function of the volume V. The chem- 
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Fig. 2.4. (a) The mean period of oscillation of the chemical system ^2.1\i - ^2.2\i as a function 
of the parameter fciij- The other parameters are given by Ig. 511 and V = 40. Each blue circle was 
obtained as an average over 10,000 periods of the stochastic model. The mean period of oscillation 
of the ODE model \2.3^ - ^2^\ is plotted as the red line, (b) The period of oscillation as a function 
of the volume V . We use k-^^ = Kj = 12.2. Parameters k2d, • ■ ■ , k^d o,re given by i2.5t . 



ical Fokker-Planck equation which corresponds to the chemical system (j2.ip ~ (|2.2p 
is a two-dimensional partial differential equation. In particular, some parts of its 
analysis are more technical than the analysis of one-dimensional problems. To get 
some insights into our approach, we analyse a one-dimensional chemical switch in the 
following section. We will present only methods which are of potential use in the 
higher-dimensional settings. A generalization of the analysis to the chemical SNIPER 
problem is shown in Section [S] 

3. One-dimensional chemical switch. We consider a chemical X in a con- 
tainer of volume V which is subject to the following set of chemical reactions (such a 
system was introduced by Schlogl [12]) 



X, 



2X 



3X. 



(3.1) 



Let X{t) be the number of molecules of the chemical X. The classical deterministic 
description of the chemical system (13. ip is given by the following mean-field ODE for 
the concentration x = X/V: 



dx 
di 



kid - k2dX + ksd X 



, ~3 



(3.2) 



To obtain the stochastic description, we first scale the rate constants with the appro- 
priate powers of the volume V, by defining 



ki = kidV, fc2 = k 



2d, 



, ksd , k^d 

ks = T^, fc4 = 



(3.3) 



V ' ■"■* y2 ■ 

Then the propensity functions of the chemical reactions (j3.ll) are given by 

Q;i(a;) = fci, a2{x) = k2X, a^ix) — k^xix — 1), a4{x) = k4,x{x ~ l){x ~ 2), (3.4) 

i.e. the probability that, given X{t) — x, the i-th reaction occurs in the time interval 
[t,t + di) is equal to ai{x) At. Given the propensity functions (|3.4p . we can simulate 
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Fig. 3.1. (a) Time evolution of X obtained by the Gillespie SSA for the chemical system l| ,?. J |l . 
The values of the rate constants are given by 113.51 1. (b) Stationary distribution of i3.lt obtained 
by the Gillespie SSA (yellow histogram) for the parameters 113.51 1. The blue curve is the stationary 
solution iS.llt of the chemical Fokker-Planck equation. 



the time evolution of the system (|3.ip by the Gillespie SSA f5] . We choose y = 1 in 
what follows, i.e. ki = kid, i = 1,2,3,4. In Figure IXTT a). we plot illustrative results 
obtained by the Gillespie SSA for the following set of rate constants 



fci ^ 2250, 



37.5, 



0.18, 



fc4 = 2.5 X 10" 



(3.5) 



We see that the system (13.11) has two favourable states for the parameter values 
p.5p . We also plot stationary distributions (yellow histograms) obtained by long 
time simulation of the Gillespie SSA in Figure lOT b). The chemical master equation 
corresponding to p.ip can be written as follows 



d_ 

m 



4 

p{x,t) =^ \a,{x + {~lY)p{x + {~l)\t) ~ ai{x)p{x,t) 



(3.6) 



where p{x,t) is the probability that X{t) — x, i.e. the probability that there are x 
molecules of the chemical species X at time t in the system. The stationary solution 
of this infinite set of ODEs can be found in a closed form TD], i.e. one can find an 
exact formula for the stationary distribution plotted in Figure [STlT bV However, our 
goal is not to solve the well-known Schlogl system using its master equation. We want 
to motivate the approach which is used later for the chemical SNIPER problem, which 
is based on the approximate description of chemical systems given by the chemical 
Fokker-Planck equation [7], see Appendix \X\ To write this equation, we have to 
consider p{x,t) as a function of the real variable x, i.e. we smoothly extend the 
function p(x,t) to non-integer values of x. Using Appendix \K[ the chemical Fokker- 
Planck equation for the chemical system p.ip is 



1^7(2^'*) = ir^{d{x)p{x,t)] - --[v{x)p{x,t) 



dt'^ ^ ' dx^ y \ /f\ J /J Q^^ 
where the drift coefficient v{x) and the diffusion coeficient d(x) are given by 



(3.7) 



,[x) =^(-l)'+^aj(x) = fci - k2X + k^xix - I) - kix{x - l){x - 2), (3. 
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1 v-^ ,s ki+k2X + k3x{x-l) + k4x{x-l){x-2) 
d[x) = -}_^a,{x) = . (3.9) 



1=1 

The stationary distribution Ps{x) = limf_+ooP(2;, t) is a solution of the stationary 
equation corresponding to (|3.7[) , namely 

da; 2 



(d(x)P,(x)) - A(„(a;)p^(^)) ^ 0. (3.10) 



Integrating over x and using the boundary conditions Ps{x) —i- as x —f ±cxd, we 
obtain 



c 
a(a;) 



^dz' 
d{z) 



(3.11) 



where c is the constant given by the normalization Jjj Ps{x)dx — 1. The function Pg 
is plotted in Figure [XlT b) for comparison as the blue line. We see that the chemical 
Fokker-Planck equation gives a very good description of the chemical system p.ip . 

The Fokker-Planck equation (|3.7p is equivalent to the Langevin equation (Ito 
stochastic differential equation) 



dX = v{X) dt + y/2d{X)dW 

where dW represents white noise. In particular, the deterministic part of the dynamics 
is given by the ODE dX/dt = v{X). Substituting dS^H]) for v{X), using (1531) and 
dividing by V, we obtain the following ODE for the concentration x = X/V: 

-7- = kid- k2dx + ksdx ix- —j - k4dX ix- — j ( ^ - 77 ] • (3-12) 

This equation slightly differs from the classical mean-field deterministic description 
p.2p , but the equations p.2p and (|3.12p are equivalent in the limit of large X in which 
X ^ ^/V. This can also be thought of as the limit of large V after a suitable non- 
dimensionalization. Let us note that the chemical Fokker-Planck equation is actually 
derived from the Kramers-Moyal expansion in the limit of large V by keeping only 
the first and the second derivatives in the expansion 7J. 

3.1. Mean switching time. Let xfi and Xf2 be favourable states of the chem- 
ical system p.ip which we define as the arguments of the maxima of Ps (given by 
p. lip ). Let Xu be the local minimum of Pg which lies between the two favourable 
states, so that Xfi < Xu < Xf2- We can find the local extrema x/i, Xu and x/2 of Pg 
as the solutions of P'g{x) = 0, which, using (j3.1ip . is equivalent to the cubic equation 

v{x)-d'{x) =0. (3.13) 

Another way to define the favourable states of the system is by considering the sta- 
tionary points of the ODE (|3.12p . i.e. the points where the drift coefficient is zero: 

v{x) ^ki- k2X + k3 x{x - 1) - fc4 x{x - l){x - 2) = 0. (3.14) 

This cubic equation has three roots which we denote yji, y^ and 2//2, where y/i < yu < 
yf2- Let us note that xji ^ y/i, Xu 7^ Vu and Xf2 7^ 2//2 because the equation (|3.13p 
differs from the equation p.l4p by the additional term d'{x). In Figure [321 a), we 
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Fig. 3.2. (a) The local extrema of the potential (i.e. the solutions Xfi < Xn < Xf2 of the 
equation i3.13l ) and the solutions yji < Vu < 2//2 "/ i'^E equation v{y) = as functions of the 
parameter ki. The critical values K and Kx are plotted as the blue dotted lines, (b) Mean exit 
times t{z) to leave the interval (— oo, (xu +3;/2)/2), given that initially X{0) = z. We compare the 
results obtained by the Gillespie SSA (yellow histogram, averages over W^ exits for each z) with the 
results obtained by the formula \3.18\ (blue curve). Parameter values are given by d,?. 511 . 



plot y/i, Hu, y/2 as functions of fci. The other parameter values are chosen as in (j3.5p . 
We also plot a:/i, a;„, Xf2- We see that there are two critical values of the parameter 
fci, namely K = 2429.35 and K^ = 2484.39. If h = K, then y/i = y„, while if 
fci = Kx, then Xfi — Xu, i.e. the number of local maxima of Ps changes at fci — K^. 
From our point of view, the more interesting critical value is K, corresponding to 
the bifurcation value of the deterministic description (j3.12p . but it is important to 
note that this value differs from the value K^. To find the value of K, we solve the 
quadratic equation v'{x) = 0. The relevant root is given by: 



Vk 



. = 1 + ^ - 77^ V(6fc4 + 2fc3)2 - 12A:4(2A:4 + h + fca). (3.15) 

0/C4 D/C4 



The value of K is determined by the equation j/ji = y„ = yk, giving 
K = kiykivk - l)iyk - 2) - kzykiyk - 1) + k2yk- 

For the parameter values (|3.5p . we obtain yk = 152.45 and K = 2429.35. 

Let T{y) be the average time to leave the interval (—00, b) given that we start 
at X{Q) ~ y. Using the backward Kolmogorov equation (adjoint equation to the 
Fokker-Planck equation (|3.7p ). one can derive a differential equation for t(j/) (see [5] 
and also Appendix [B|) . namely 



) — 
dy 



dV 



v{y)—{y) + d{y)^ 2 



(y) for 2/ e (-00,6), 



(3.16) 



with boundary conditions 



dr 
Ay 



(-00) = 0, r(&) = 0. 



(3.17) 
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Solving (|3.16[) . we obtain 

— (y) = -exp 
ay 



v{z) 



1 



diz) 

V 



d{x) 



exp 



w(z) 
diz) 



dx 



d{y)Ps{y) J- 



Psix)dx, 



where Ps is the stationary distribution (|3.1ip . Integrating over y in the interval (z, b) 
and using (|3.17p . we obtain 



t{z) 



dr 
dy 



{y)<iy 



d{y)Ps{y) y_ 



Ps{x)dxdy. 



(3.18) 



The important characteristic of the system is the mean time T(ki) for the system to 
switch from the favourable state Xfi to the favourable state Xf2- To determine this 
we set 



b = 



_ f{x.,+Xf2)/2, 

284, 



for ki < Ka:, 
for ki> K^, 



(3.19) 



and we define T[ki) precisely as the mean time to leave the interval (— oo, b) given that 
X(Q) = 0. This definition provides a natural extension of the mean switching time 
for values fci > Kx- Note that the number 284 in (|3.19p is the value of {xu + a;/2)/2 
at the critical point ki — Kx- Using (|3.18p . we can compute T{ki) as the integral 



r(fci 



exp 



v{z)_ 
d{z) 



Az 



1 



-oo d{x) 



exp 



v(z)_ 
d{z) 



dz 



dxAy 



(3.20) 



where b is given by p.l9p . Let us note that we could also put the right boundary b at 
the points x„ or y„. If the number of molecules reaches the value a;„ (resp. yu)^ then 
there is, roughly speaking, 50% chance to return back to the favourable state xji and 
50% chance to continue to the second favourable state Xf2 ■ Thus the mean switching 
time between the states could be estimated by multiplying the time to reach the point 
Xu (resp. yu) by the factor of 2. The problem with this definition is that Xu ^ yu- 
The difference between the results obtained by choosing the escape boundary at x„ or 
yu is not negligible because the drift is small between Xu and yu- On the other hand, 
the drift is large at {xu + a;/2)/2. Replacing [x^ + a^/2)/2 by any number close to it, 
e.g. by {yu + j//2)/2, leads to negligible errors. In that sense, the boundary b given by 
p.l9p yields a more robust definition of the mean switching time. The other reason 
for presenting analysis for the definition (|3.20p is that it is naturally transferable to 
the chemical SNIPER problem which is the main interest of this paper. 

In Figure [3?2r b). we compare the results obtained by the Gillespie SSA and by 
formula (|3.18p for the parameter values (|3.5|) and the right boundary at 6 = {xu + 
a;/2)/2 = 310.8. We have Xfi = 99.9, Xu = 225.6 and x/2 = 396.0 for the parameter 
values p.5p . In Figure [3?^ b). we plot mean exit times computed as averages over 10^ 
exits from the domain [0, 6], starting at X{Q) — z for every integer value of z < 6. The 
results are in excellent agreement with the results computed by the formula (|3.18p . 

In Figure [3?2r b). we also see that there is no significant difference between the 
exit times if we start at X{Q) = Xfi or X{0) = for fci ^ Kx- This justifies the 
choice X{0) = in the definition of the mean switching time T(fci). If fci w Kx, than 
Xfi is close to Xu and the results obtained by the starting point at X(0) = and at 
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X{Q) — Xfi will differ. In the definition (|3.20p . we use X{0) — for any value of fci 
because (i) the results for X{Q) = will provide some insights into the estimation 
of the period of oscillation of the SNIPER problem studied later; (ii) the results are 
robust to small changes in the initial condition X{0) — 0; (iii) the starting point 
X{Q) — is defined for any value of fci (note that Xfi is not defined for fci > K^)- 

In the remainder of this section, we provide approximations of T{ki) defined by 
the formula (|3.20p . We fix the parameters ^2, ^3 and ^4 as in ()3.5|) and we vary 
the parameter fci. In Section [321 we provide the estimation of T{ki) for ki ^ K 
(outer solution). In Section [3. 31 we provide the estimation of T(fci) for ki~K (inner 
solution) . Finally, in Section 13. 4[ we match the inner and outer solutions to obtain 
the uniform approximation for any fci. 

We note here that one could approximate T(fci) by approximating the integral 
(|3.20p . Using the method of steepest descent, we would obtain the generalization of 
the well-known Kramers formula 

. s . X 27rexpf$(a;„) — <J)(a;fi)l 
T fci ~ Tfc fci = \ ^ ' ^Jlll- 3.21 

where ^{x) is defined by $(a;) — — f^ v{y)/d{y)dy + log{d{x)). This approximation is 
valid for ki <^ K |TT]. Haataja et al. [S] suggest another generalization of the Kramers 
formula replacing (i(x„) in the denominator of (|3.2ip by {d{xu) + '^(2;/i)}/2, but this 
approximation is worse than (|3.2ip . Unfortunately, such an integral representation is 
not available for higher-dimensional Fokker-Planck equations, including the SNIPER 
problem studied later. Since our main goal is to present methods which are applicable 
in higher-dimensions also, in Sections I3.2[ 13.3] and 13.41 we focus on approximating the 
mean switching time T(ki) by analysing the r-equation (I3.16p . The chemical Fokker- 
Planck equation and the corresponding r-equation are available in any dimension - 
see Appendix |A] and Appendix |b1 

3.2. Approximation of the mean switching time T(fci) for fci <C K. For- 
mula (j3.9p and the parameter values (13. Sp imply that d{x) is of the order 10^ while x 
is of the order 10^. Considering small e ^ 10"^, we use the following scaling 

V , d y /„ „„x 

^ = ^^ ^=^' y^-,' ^ = ^^- (3-22) 

Then we can rewrite the equation (|3.16p to 

v{y)^iy) + ediy)^iy) = -i- (3.23) 

Later we will see that in fact r is exponentially large in e, so that the left-hand side of 
this equation dominates the right-hand side. Then r will be approximately constant 
except when w « 0, that is, the main variation in r occurs near y = y^- For y close 
to y^, we use the transformation of variables y — yu + \/^'7 and approximate 

d^ 

v{y)~V~ev—{y^), d{y)^d{y^), r»l, (3.24) 



dy 



to give 



dw , , dr , , — , , d T 



^d^^^-^d^^^^ + ^^^-W^''^"'^" ^^-^^^ 
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Integrating over ry, we get 

exp 



1 diJ 



(fJ^' 



Cl 



dr 



d'/ d(2/.J 



(3.26) 



where Ci is a real constant. To determine ci we cannot use the local expansion near 
y^ alone; Ci is determined by the number —1 on the right-hand side of p.23p . which 
is the only place the scale of r is set. To determine Ci we use the projection method 
of Ward [2] . The stationary Fokker-Planck equation (|3.10p can be rewritten 



- ^ [vrnPsin + e^ \d{v)PM] = 0, 



(3.27) 



which is the adjoint equation to the homogeneous version of (|3.23p . Multiplying the 
equation p.23p by Ps, the equation (|3.27p by r, integrating over y in the interval 
[—oo^y^ and subtracting the resulting equations, we obtain 



-T^[Ps{y)v{y)f{y)]dy 
-oo dy 



+ £ 



Vu 



_d_ 
dy 



Ps{y)d{y)'^{y) 
dy 



_d_ 
dy 



1y)^ [d{y)Psiy)] 



dy = - 



Ps{y)dy. 



Evaluating the integrals on the left hand side with the help of (|3.27p and the boundary 
condition p.l7p . we find 



dy 



1 



ePsiyJdiyu) 



Ps{y)dy. 



(3.28) 



Since Ps is exponentially localized near y^^ and y^, Ps{yu) is exponentially small in 
s. Using (p:^ in dS^S]), we find 



Cl =%J^(0) = V-sd{yu)^{yJ - ~^^^) 



so that 



dr 
drj 



1 



exp 



/iP,(yJd(yJ 
Integrating over rj in [— cx),oo], we obtain 

V2^ 



lim T{'if) 

rj — ^ — oo 



■Ayu)Ps{yu) 



1 diJ 
"2rf(yjd^ 



dv 



(yuW 



■1/2 



PMdy, 



Ps{y)dy- 



Ps{y)dy, 



(3.29) 



where we used that t(?7) -^- as 77 — » 00, i.e. the switching time is going to zero if 
the starting point approaches the boundary b, see Figure [32Ib). The limit rj —^ ~oo 
on the left hand side of (I3.29P is an approximation of the mean switching time T(fci). 
It is the plateau value of r(z) in Figure [HT^ bV Transforming (|3.29p to the original 
variables, we obtain the following approximation for the mean switching time 



r(fci)«ra(fci) 



2tt 



y/d{yu)Psiyu) Vdy 



dv 



-{Vn 



-1/2 



Ps{y)dy. 



(3.30) 
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(a) 



2150 



_T(k ) exact 

T (k,) outer 
T, (k ) Kramers 




2200 



2250 



2300 



(b) 



T(k ) exact 

T (k ) outer 

- 0^ 1' 

I Kramers 




2500 



Fig. 3.3. Approximations H3.21\i . H3.30\i and l|3.3J|l of the mean switching time T{k\). The red 
line is the exact value given by ^3.20\ . (a) Region k\ <^ K for which the derivation was made, (b) 
Behaviour of approximations close to the critical points ki = K and fci = Kx- The critical values 
K and K^ are plotted as the dotted lines, K < Kx ■ 



Let US define the potential ^ by 

*(a;) 
Then p.30p can be rewritten as 



d{z) 



Ta(fci) = ^/2TTd{yu) exp [*(?/„)] ( -r-iVu 



du , 



dz. 



-1/2 



exp[-^(^)] 
d{v) 



Ay. 



The local extrema of the potential 4" are equal to y/i, ?/„ and y/2- Using the Taylor 
expansion ^(y) « ^(?//i) + {v ~ V fiY^" {v fi) /"^ in the integral on the right hand side, 
we approximate 



Ay 

Tih) « ra(fci) « eM^(.yu)-^{yfi)] [ ^ivu) 



-1/2 



2TTy^d{yu) 



d{yn).JW^, 



Using the definition of ^, we obtain the following approximation of the mean switching 
time 



T(fci) « To{ki) = 2Try/d{yu) exp 



v{z) 



dz 



dv , 
dy 



iVu) 



dv 
dy 



iVfi) 



diVfi) 



-1/2 



(3.31) 
We will call To{ki) the outer solution. In Figure [5?5f a). we compare the approxima- 
tions (|3.30p and (|3.3ip with the exact value given by (|3.20p . We see that T^ and To 
provide good approximation of the mean switching time. In Figure I^TST b). we present 
the behaviour of approximations close to the critical point fci — K . They both blow 
up at the point fci = K. We also plot the results obtained by the Kramers approxima- 
tion Tfe(fci) given by (I3.2ip . We again confirm that it provides a good approximation 
for fci ^ K, see Figure lXST a). but it blows up at the point ki = Kx, see Figure [373l b). 
In the next section, we present the inner approximation which is valid close to the 
critical point fci = K. 
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3.3. Approximation of the mean switching time T(fci) for ki w K. If 
ki = K, then y^ = yt^ = y^. where y^ is given by p.lSp . Consider the drift coefficient 
given as a function of y and the parameter fci , namely 

v{y,ki) = fci - k2y + k3y{y- 1) - k^yiy - l){y-2). (3.32) 

If y is close to y^, and ki is close to K , we can use the Taylor expansion to approximate 



v{y,ki 



2 



/_ _ n2 



;;^{y,,K){y-y,r + {k,-K), 



(3.33) 



where we used 

Qv _ dv _ d V _ d V _ 

— iy„K) = l and v{y,,K)^-{y„K)^^=^^{y,,K)^^^{y,,K)^0. 

We use the transformation of variables y = yj. + e^^'^ rj and ki = K + e^l'^ k. Then 
(|3.33p implies 



tJ(y, fci) « 6^/3 ( \-l{y^^K)rf + K 



(3.34) 



Now in p.22p we need to replace t — er hy t ^ e'^/^ t. Then ()3.23p reads as follows 



-l/3^fc' 



dr 



=-2/3jfc- 



d?T 



v{y)^{y)+e'''d{y)^iy) = -l. 



(3.35) 



Using (I3.34P and approximating d{y) ss d{yf.), we get 






-1. 



(3.36) 



This time, since we are close to the bifurcation point, the right-hand side is not 
negligible. Equation p.36p can be rewritten as 



— (^d(y,)exp 



i.(fc(w^<"-?-" 



dr 
drj 



' exp 



1 /I d^v 77^ 



Integrating over 77, we get 
dr 1 



d^ c?(yfe) 



exp 



1 fld^v,_ .rf 



exp 



1 (id^v e ^ , 



de 



(3.37) 



We chose the exiting boundary so that T{r]) — > at ry ^ 00. Consequently, integrating 
p.37p over 77 in [—00, 00], we obtain 



Km r(?7) = =— - / / exp 



1 



rf(yfc) V% 



q2 rt' — F^ 



6 



d^drj. 
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Substituting 

V3d(2/fe) V 
we obtain 

where 



-1/3 



[u + v], i = 



lim t(77) 



TT \fl2 



rf(^J/3^ 



1 d^v,_ ^, 



7J(/3k) . 



-1/3 



and the function H 



/3 = -24V3rf(yJ-V3(^g(y^,^)y 
M ^ (0, oo) is given by 

H{z)= ^^ ^ ^-dv. 

Jo Vv 



-1/3 



(u-v), 



(3.38) 



(3.39) 



(3.40) 



The limit of ^(ry) as 77 — > —00 is the approximation of the mean switching time T(ki) 
for fci close to K. Transforming back to the original variables and using (|3.39p . we 
obtain the inner solution Ti[ki) in the following form 



dy 



2iyk,K) 



xi/(-[fci-if]24i/3(%fe))-2/3 



-2/3 



dy 



,{yk,K) 



-1/3N 



(3.41) 



In Figure [5^ a). we compare the approximation (|3.4ip with the exact value T[ki) 
given by (|3.20p . We also plot T{xfi) given by (|3.18p for b = (a:„+a:/2)/2- As discussed 
before, t(x/i) could be considered as another possible definition of the mean switching 
time. We see that the value of the inner approximation Ti{ki) at the critical point 
ki — K lies between the exact value T(fci) and T{xfi). Thus we confirm that Ti{ki) 
is a good approximation of T(fci) for fci close the the critical point K. 

3.4. Uniform approximation of the mean switching time Tiki). In Sec- 
tion l3.21 we obtained the outer approximation (I3.3ip of the mean switching time which 
is valid for fci <C K. In Section [5T51 we obtained the inner approximation (|3.4ip of the 
mean switching time which is valid for fci w K. In this section, we match the inner 
and the outer solutions to derive the uniform approximation of the mean switching 
time. First, we investigate the asymptotic behaviour of H{z) as z — > 00. Let z > 0. 
Using the substitution w = z^^l'^v in the definition (|3.40p . we obtain 



F(z) 



= zl/4 



1 



exp 



z3/2 U 



Aw 



n 
— exp 



/o Vw 
Consequently, the outer limit of the inner solution p.4ip is 

v/|fci-X| [dy 



2^3/2' 



T- — 



TT^ly'^'^) 



exp 



4^/2 
3d(2/fe) 



3\/3 



92. 



as z 



fci-^r(i^(^^^^) 



-1/2- 



(3.42) 
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(a) 



(b) 
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:e 5 

o 
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S3 







T(k^) exact 






_T^(k^) outer 






T.(k ) inner 






-VS) 


\\\ \ 


.-x(x„) 


"^ ^ ^^v^^*^^"— ^^^' 


- 


~^^^5v. 


^^^trr— -^^as-^ 





2300 2350 2400 2450 2500 2550 2600 




2600 



Fig. 3.4. (a) The inner solution \3.41\ as a function of k\ (black line). The red line is the exact 
value T(k\) given by 1 13.201 1. We also plot the outer solution 1 3. 3 It (blue line), the outer limit of the 
inner solution ]3.4S[ (black dashed line) and T{xf-i) given by i3.18l for b = (xu +Xj^2)/2 (red dashed 
line), (b) The uniform approximation T^nifi^i) which is given by \3.43\ (b lue line). We also plot 
the exact value T(k\) given by ^3.201 (red line) and T{xfi) given by i3.18[ for b = {xu + Xf2)/'2 
(red dashed line). 



In Figure [3^ a). we plot Ti-oiki) as the black dashed line. It is easy to check that the 
inner limit of the outer solution To(fci) is equal to the outer limit of the inner solution 
Ti{ki), i.e. it is also given by Tj.o(fci). Thus we can define the continuous uniform 
approximation of the mean switching time T{ki) by the formula 



-Lunifykl) 



^ I Toih) + T,{ki) ~ r,;;o(fci) - C, for fci < K- 
~lT,(fci), forfci>if; 



(3.43) 



where the constant C is defined by 



C 



lim {To{ki) ~ T,,o{ki)) . 



Thus, for fci < K, we add the inner and the outer solutions and subtract the "overlap" 
solution (the inner limit of the outer solution) which has been double-counted. In order 
to make the approximation continuous at ki — K , we also subtract the constant C 
(which is in fact a higher order term). We approximate T{ki) by the inner solution 
for fci > K. In Figure [5^ b). we plot the uniform approximation Tunif{ki) together 
with the exact value Tiki) given by (|3.20p . We also plot T{xfi) given by (|3.18[) for 
h = (xu + a;/2)/2. The comparison is excellent. 

4. Numerical results obtained by solving the chemical Fokker-Planck 
equation for the chemical SNIPER problem. We consider the chemical system 
(|2.ip - (|2.2p . Substituting the propensity functions (|2.7p into the equation (jA.ip in 
Appendix [XI we obtain the chemical Fokker-Planck equation 



dP 

'dt 






[^^^l + ^C'^-^l + ^t^^^] 



d 



d 



,,K^]-^K^] 



(4.1) 
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where P{x,y,t) is the joint probabiUty distribution that X{t) = x and Y{t) — y, 
X e [0, cx)), y g [0, cx)), and the drift and diffusion coefficients are given by 

Vx{x, y) = fc2j/ - k^x{x - l){x - 2) + kix{x - 1) - fcsx, (4.2) 

Vy{x,y) = -kYx{x - l)y + k^xy - k2y + ki, (4.3) 

dx{x, y) = [k^y + k^x{x - l)(a; - 2) + kix{x - 1) + fc3x]/2, (4.4) 

dy{x, y) = [kjx{x - l)y + k^xy + k2y + fci]/2, (4.5) 

dxy{x,y) = -k2y. (4.6) 

The stationary distribution Ps{x, y) = hm^^oo Pix, y, t) can be obtained as a solution 
of the corresponding eUiptic problem: 

- £ [^^ ^^] + &y [^- ^^] + W f'^ ""'^ - i f^^ ""'^ - I f^^ "^^^ ' ^'-'^ 



f-OO /'OO 

/ / -Ps(a;,y)dxd2/= 1, 
Jo Jo 



(4.8) 



Ps{x, y) > 0, for all (x, y) G [0, oo) x [0, (»). (4.9) 

To solve the problem (|4.7p - (|4.9p approximately by the finite element method (FEM) 
we have to reformulate it as a problem in a finite domain. Therefore, we first restrict 
the domain [0, oo) x [0, oo) to a rectangle S which has to be sufficiently large to contain 
most of the trajectories. For example, Figure [2?2r b) shows the illustrative trajectory 
for the parameter values ()2.5|) and V — 40. In this case, the rectangle S can be chosen 
as S* = (0, 500) X (0, 2000). On the boundary dS we prescribe homogeneous Neumann 
boundary conditions. We reformulate (|4.7p - (|4.9p as the following Neumann problem 
in the domain S 

- div(y^VP, + Psh) =0 in S', (4.10) 

{AVPs + Psh) • n = on as*, 

where V stands for the gradient, n denotes the outward unit normal vector to dS, 
and the 2x2 symmetric positive definite matrix A and the vector b are given by 



_ , -dx -dxy/2\ y^_ f _ 9dx _ 1 ddxy _ ddy ^ 1 dd^y 



T 

^-dxyll -dy }' ""' \^^'^'2~d^'^'"~^''2~d^ 

(4.11) 
To define the FEM solution, we need the weak formulation of ()4.10p . The weak 
solution Pg G H^{S) is uniquely determined by the requirement 

a{P,,ip)=0 V^eif^(5), 

where H^{S) stands for the Sobolev space W^'^{S) and the bilinear form a(-,-) is 
given by 



a{P,,^) = / {AS/P, + Psh) ■ V^dxdy. 
Js 
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Solution of Fokker-Planck Equation 



Stochastic Simulation 




Fig. 4.1. (a) The FEM solution of the Fokker-Planck equation OTTOl restricted to the subdomain 
(0,80) X (550, 1150). (b) Ps{x,y) obtained by the Gillespie SSA. We use the parameter values 112.51 1 
and V = 40. 



We use first-order triangular elements. First we define a triangulation T^ of the 
domain S and the corresponding finite element subspace Wh of the piecewise linear 
functions: 



Wh = Wh e H\S) : ^„|k £ PHk), K g %}, 



(4.12) 



where P^{K) stands for the three-dimensional space of linear functions on the triangle 
K € Th- The finite element problem then reads: find Pg.h G Wh such that 



aiPs,h, Vh) = ^iph G Wh- 



(4.13) 



Both problems (|4.10p and (|4.13p posses trivial solutions Pg = and Pg.h = 0. To get 
a non-zero solution we use appropriate numerical methods for finding eigenvectors 
corresponding to the zero eigenvalue. The computed nontrivial solution Ps^h is then 
normalized to have J^ Ps,h{x,y)dxdy = 1. In Figure l4?lT a). we present the FEM 
solution Pg^h to the problem (|4.13|) obtained on a uniform mesh with 2^® triangular 
elements. We use the parameter values (|2.5|) and V = 40. We plot only the stationary 
distribution in the subdomain [0, 80] x [550, 1150], i.e. the part of the X-Y phase space 
where the system spends most of time. Running long time stochastic simulations, 
we can find the stationary distribution by the Gillespie SSA, which is plotted in 
Figure l4?lT b) for the same parameter values (|2.5p and V — 40. The comparison with 
the results obtained by the chemical Fokker-Planck equation is visually excellent. 
Plotting the numerical results in the whole computational domain S* = (0, 500) x 
(0,2000), we do not see any additional information because most of Pg is localized in 
the subdomain shown in Figure [JT] We have to plot log(Ps) to see the underlying 
"volcano-shaped" probability distribution - see Figure [4?2] (a). Let us note that there 
is no bifurcation in the features of Pg, i.e. the probability distribution is "volcano- 
shaped" both before and after the ODE critical value. 

4.1. Computation of the period of oscillation. Let us consider the SNIPER 
chemical system (|2.1|) - p.2|) with the parameter values (|2.5|) and V — 40. An illus- 
trative stochastic trajectory is shown in Figure [5T21 Let the domain Q. be defined by 
^ = {[x,y\\ X < 200}. In each cycle (as we defined it on pageO, the trajectory leaves 
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(a) (b) 
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X 

Fig. 4.2. (a) Logarithm of the FEM solution of the Fokker- Planck equation OTTOl in the 
domain S = (0,500) X (0,2000). (b) The solution T(x,y) of H4.15^ computed by the FEM. We use 
the parameter values 1 12.511 and V = 40. 

the domain Q. However, the time which the trajectory spends outside the domain 57 
is much smaller than the time which the trajectory spends inside fl. This observation 
is confirmed by Figure [?^ a). The probability that the system state is outside J7 is 
several orders of magnitude smaller than the probability that it lies inside fi. Thus 
the period of oscillation can be estimated as the mean time to leave the domain fi 
provided that the trajectory "has just entered it". Let t{x, y) be the average time to 
leave the domain fi provided that the trajectory starts at X{0) ~ x and Y{0) — y. 
As shown in Appendix [Bl T[x,y) evolves according to the equation (jB.SP : 

'^^9^+''-9^ + '^9^+"^9^ + "^9^ = -'' f-[--'^]^"' (4.14) 

together with the boundary condition r(200, y) = for y g R. To solve this problem 
approximately, we truncate the domain Q to get the finite domain S = (0, 200) x 
(0,2000). We consider homogeneous Neumann boundary conditions on the parts of 
the boundary which are not on the line x — 200 and we rewrite the problem (|4.14p in 
the form 

-div(y^Vr)+b- Vr = -1 in S^, (4.15) 

T = on the line x = 200, 
(AVt) • n = on the lines y = 0, y = 2000, a; = 0, 

where the 2x2 matrix A and the vector b are given by ()4.1ip and the coefficients 
dx, dxy, dy, Vx, and Vy are given by (I4.2|) - (I4.6|) . Notice the difference between ()4.15p 
and (|4.10|) and that S* C 5". As in the previous section, we define the weak solution as 
T ^W satisfying 

a{T, (p) — / — 1 • ^<^xdy ^Lp G W ^ 
Js 

where W = {v e H^{S) : v = Q on the line x = 200} and 

a{T^ ip) = / AS/t ■ \7ip dxdy + / b • \7Tip dxdy. 
Js Js 
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The FEM solution r^ G Wfi is defined by the requirement 

a{Th, fh) = / -1 • fh dxdy, Viph € Wh, 

Js 

where Wh C M^ is the space of continuous and piecewise linear functions based on a 
suitable triangulation T^ of 5" (compare with (|4.12[) ). The triangulation used does not 
follow the anisotropy of S and it consits of elements close to the equilateral triangle. 
It has about 500, 000 triangles and the number of degrees of freedom is about 250, 000. 
The FEM solution r^ is shown in Figure HTW b) . Using the computed profile, we can 
estimate the period of oscillation. One possibility is to find the maximum of the sta- 
tionary distribution Ps(x,y) which is attained at the point [x/,?//] = [29.30,781.25], 
and estimate the period of oscillation as Th{xf,yf) — r^ (29.30, 781.25) = 134.1. This 
compares well with the period of oscillation estimated by long stochastic simulation 
for the same parameter values. We obtained 130.4 as an average over 10^ periods. 

The estimation of the period of oscillation is analogous to the estimation of the 
mean switching time of the Schlogl problem (|3.ip which was studied in Section [3] We 
introduced the boundary b given by (|3.19p and we asked what the mean time to leave 
the domain (— oo, b) is. The biggest contribution to the mean exit time was given 
by the behaviour close to the point Xu- The SNIPER equivalent of the boundary b 
(resp. point Xu) is the line x = 200 (resp. the saddle). In the case of the Schlogl 
problem (|3.ip . we saw that the difference between the escape times from X{0) — 
and ^(0) = Xfi increased as we approached the bifurcation value. Similarly, the 
estimation of the period of oscillation of the SNIPER problem by t/i (a;/ , y/ ) is fine for 
kiii <C Kii, but it might provide less accurate results if the saddle and the stable node 
of the SNIPER problem are too close, i.e. if /ci^j is close to the bifurcation value K^i. 
Thus it is better to estimate the mean period of oscillation as the mean time to leave 
the domain fl, starting from the subdomain j d Q where stochastic trajectories enter 
the domain fl. We approximate the mean period of oscillation as a weighted average 
of r(a;, y) over a suitable subdomain 7, namely by 



T{x,y)Ps{x,y)dxdy 

ni) = '-^^ . (4.16) 

Ps{x,y)dxdy 



7 

Since the trajectories are entering the domain S at the smaller values of Y and leaving 
the domain at the larger values of y, it is reasonable to choose the subdomain 7 as 
the line segment parallel with the Y axis between the X axis and a suitable treshold 
value. Three choices of 7 are shown in Figure I4.3f a) as thick black, magenta and 
green lines. The black line is the line from y = to the unstable node. The magenta 
line is a segment oi x = 130 and the thick green line is a segment of a: = 100. In the 
case oi X — 100 (resp. x = 130), we define the threshold value as the intersection of 
the line x = 100 (resp. x — 130) with the stable direction at the saddle point. Note 
that we use the nuUclines and the saddle point of the system of ODEs 



dx dy 

— — = Vx(x,y), — - 

di ^ '^^' dt 



= vx{x,y), —=Vy{x,y) (4.17) 



which slightly differs from the classical mean-field ODE description (I2.3p - (|2.4p . The 
ODE system (HTTl) is equivalent to the ODE system ([^ - ((^ in the limit F ^ 00. 
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Fig. 4.3. (a) NuUclines and steady states of the ODE system | |^.i7| l together with the lines 7 
which are used in the formula (RTTfill. The lines 7 are plotted as the thick green line (for x = lOOj, 
the magenta line (for x = 130 j and the black line (for the line to the unstable node), (b) The period 
of oscillation as a function of the parameter ki^ . The results obtained by the Gillespie SSA (blue 
circles) are compared with the results obtained by (pTTCk for the lines 7 shown in the panel (a). 



See also the discussion after the equation (|3.12p about the differences between p.2p 
and (Pl^ . 

In Figure H3f b). we present estimates of the period of osciUation computed by 
(|4.16p as a function of the parameter kid- The results obtained by the Gillespie 
SSA (blue circles) and by the ODEs (|2.3p ~ (|2.4p (red line) were already presented 
in Figure I^TW a). The green, magenta and black curves correspond to the results 
computed by (|4.16p for the lines 7 of the same colour in Figure[4?3Ka). We also present 
results computed by T{xf,yf) where [xf,yf] is the point where Ps{x,y) achieves its 
maximum. The dotted line denotes the bifurcation values of the ODE systems (12.3^ - 
(|2.4p and (|4.17p . These values differ for finite values of the volume V. In Figure ll^ a). 
we show the dependence of the bifurcation value of kid on the volume V for both ODE 
systems. The ODE system (|2.3p - (|2.4p is independent of V, and so its bifurcation value. 
Fixing the value of kid at the bifurcation value Kd = 12.2 of the ODE model (|2.3p - 
(p:i)) . the period of oscillation of the ODE system ([^ - ([^ is infinity. The ODE 
system (|4.17l) does not have a limit cycle at all as Kd = 12.2 is below its bifurcation 
value which we denote K. However, we saw in FigureH^b) that the stochastic system 
has oscillatory solutions with a finite period for kid = Kd- The period of oscillation 
as a function of the volume V can be also computed by (|4.16p . The results are shown 
in Figure B^ b). We use the same lines 7 as in Figure [THT a). We have to take into 
account that the system volume V changes and that the phase plane axes (the number 
of species particles) scale linearly with V. Thus we define the lines 7 as segments of 
X = 0.67x", X — 0.87a;" and x — x"^ where x" is the x-coordinate of the unstable 
node. This definition gives for V — 40 the lines 7 plotted in Figure [4?3ta). We also 
have to scale the domain S with V: we use SV/AO instead of S. Fortunately, we can 
simply rescale the triangulation. Thus the volume dependence of the computational 
domain does not causes any additional problems. Let us note that the y threshold 
for the lines at x = 0.67a;" and x — 0.872;" is defined as an intersection with the 
stable direction at the saddle point of (|4.17p . As shown in Figure I4.4r a') , the saddle 
point of the ODE system (|4.17p is always defined for kid = Kd because its (volume 
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Fig. 4.4. (a) The dependence of the bifurcation value of the parameter k^^ on the volume V for 
the ODE system \-i.n\ (solid line) and for the ODE system ^2.3^ - ^274^ (dashed line). The other 
parameter values are given according to \2.5\ . (b) The period of oscillation as a function of the 
volume V . The results obtained by the Gillespie SSA (blue circles) are compared with the results 
obtained by CTTfil for the lines 7 (scaled by V) which are shown in Fiaure W^Ti a), 



dependent) bifurcation value K satisfies K > Kd- Let us note that for fci greater than 
the bifurcation value of (|4.17p . we use the lines 7 computed at the bifurcation point 
(since there is no saddle defined for such values of kid). This definition was used in 
Figure SSlJb) for fcid > K. 

5. Analytical results. In this section, we derive analytical formulae for the 
dependence of the period of oscillation of the chemical SNIPER problem ()2.ip - (l2.2|) 
on the parameter kid and on the volume V , i.e. we obtain formulae for the behaviour 
shown in Figure [FM In particular, we generalize the one-dimensional approach from 
Section [3731 to the chemical SNIPER problem. As we mentioned before, we denote the 
bifurcation value of the parameter kid of the ODE system (I4.17P by K. If kid — K, 
then the saddle and the stable node of the ODE system (|4.17p coincide at one point 
which we denote as [xk , yk] ■ Let us define 



A 




{xk,yk,K), 



where the dependence on parameters k2, . ■ ■ ,k-^ is not indicated, because their values 
are fixed and given by (|2.5p . The eigenvalues of A are Ai = and A2 = — Aq < 0. We 
denote the corresponding eigenvectors as Ui = (un, M21) (for the eignevalue zero) and 
U2 = {ui2, "22 ) (for the eigenvalue — Aq). The eigenvector Ui is in the direction of the 
limit cycle. There is an obvious separation of time scales. The dynamics is much slower 
in the direction of Ui than in the direction of U2. Thus we first transform the variables 
from the Cartesian coordinate system x and y to the directions of eigenvectors of A 
and then we analyse the transformed r-equation. To do it systematically, we use the 
following scaling (compare with (13.2211 ) 



X 

e 



-2 ' 
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Then the r-equation (|4.14p reads as follows 



_ d'^T — d'^T — d'^T _ dr _ dr 



£ d^ 7^ + £ d,y -=T^ + e dy —^ + V, — + vy — = -e. (5.1) 



For kid close to K and [x,y] close to [xk,yk], we use the local variables 77, ^ and k 
which are related to x and y as follows 



;)^©-(s: £)(-;"?)■ ^...^^^.^''■-^-»- (-> 



Then the line ^ = corresponds to the direction of the limit cycle and ?/ = to 
the stable direction. Let us denote q = [x/c , y j, , -R'] to shorten the following formulae. 
Using the Taylor expansion at the point q — \xk , y^ 1 K], we approximate 



rf.(^,y,fci)«d.(g)+ei/3ci(d,)77 + ei/2c2(4)e + £'/'c3(d,)r7' + o(e^/'), 
d^y{x,y, ki) « d^y{q) + e^/^ci(d^y) 77 + e^/^C2(da,y) £, + e^^^c3(d.^y) rf + o{e^/^), 

dyix,y, fcl) « dyCq) + e^/^Cl{dy) T] + E^ ' ^ C2{dy) ^ + E^^^'c^idy) if + o{e'>/^), 



where Ci{dx)^ Ci{dxy) and Ci{dy) are constants. To systematically derive the formulae 
for the mean period of oscillation, we need to take all the terms above into account. 
However, we will show that the results are actually independent of some coefficients 
in the expansion. To save space, we explicitly specify only those coefficients which 
actually appears in the final formulae. They are given in Appendix [C] Using the 
fact that Ui = (uii,M2i) (resp. U2 — (ui2,M22)) is an eigenvector of the matrix A 
corresponding to the eigenvalue Ai = (resp. A2 = — Ag < 0), we obtain the expansion 
of Vx and Vy in the form 



Vx{x,y, fci) « -e^/^XoUi2 C + e^/^ci(U^) i]^ + e'^^'^c^ivx) V^ 

Vy(x,y, fci) w -£^/^AoU22 C + £^^^(ci(%) 77^ + k) + e'^/'^£,C2{vy) r? 
+ e{c3{vy) e + C4{vy)) v') + e'/'c,{vy) v'^ + o{e^^'), 



where the constants Ci{vx) and Ci{vy) are given in Appendix [Cl Using (|5.2p . the 
derivatives transform as follows 



d _ d d-q ^ d dS, _ _i/3 U22 d ^_^i2 U21 d 



dx drj dx d£, dx det U drj det U d^ 

d_ _ d_&n _^^ _ _ -1/3 "12 d „i/2 "11 9 
dy di-j dy d^ dy det U dr] det U d£, ' 
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where detU — U11U22 — Ui2W2i- Substituting our approximations to (|5.ip and using 
the scaling r = e^/'^r, we obtain 



,dT 



I ^ ^9t 9^r\ i/fi / 2 



9r d^T \ 

g4 ^ ^ 1 



|2= 



a^T 












_2/3 






a- 



01^ 



crigrj 



..dr 

drj 



0-20 f? ^ 









,5r 



(5.3) 



where the coefficients ct; are given in the Appendix [C] We assume that r — > as 
rj —t 00 and that r is bounded as 77 ^ — cx) and ^ ^ ±00. We expand r as 

T ^ TO + ei/Vi + £1/^2 + ei/Va + £2/^4 + o{e^"^). 

Substituting into (|5.3p . we obtain the following equation, at 0{e^), 






(5.4) 



Integrating over ^, we obtain 






Civ) exp 



2(71 



where C is a function of 77. Since r is assumed to be bounded as ^ — > ±00, we must 
have C{ri) = 0. Consequently, we obtain 



To = To{ri), 



(5.5) 



i.e. To is not a function of ^. Using (|5.5p in (|5.3p . we obtain the following equation, 
at 0(ei/6), 



Ao C^^ - cTi 



d£, " ' ae 



= 0, 



which similarly implies 

n = Ti{r]). 

Using (|?3)) and (H^D in ((0|l . we have, at ©(e^/^), 

, .5x2 9V2 2 5^0 9to 

Ao ? ^F ~ ^1 ^^ = erg 77 ^r- + erg K-— 



CTg 



Thus 



-<Ji-^[cxp 



xoe 



2o-i 



9t2 
9^ 



exp 



xoe 



2ai 






<9ro 



(5.6) 



(5.7) 



(5.8) 
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Integrating over ^ in (—00,00), we obtain the solvability condition 



[<^5V 



aen 



dr] 



CT8 



dif 



-1. 



(5.9) 



This equation is the SNIPER equivalent of equation p.36p for the chemical switch, 
and can be solved to find the leading-order period of oscillation. Before we do so, we 
proceed to higher-order with the expansion to determine the correction term to tq. 
Using ((0|) in (|5?8l) gives 



dT2 d^T2 



0. 



Thus, again. 

Using dSH), dEll) and ((5T0| in ([Ell), we have, at 0{e^/'^) 



Ao^ 



Sra 



(Tl- 



9V3 



5C ' "' de 

which is equivalent to 
d 



CbV 



,dTi 



(JQ K- 



drj 



CT8 






cfwiri' 



drj 



CTi^ I exp 



^oe 



2cri 






= exp 



2cri 



,CT5 7? 



0-6 K I 



ari 



CTs- 



aVi 



crio C'/ 



drj ' ° 977^ 
Integrating over ^ in (—00,00), we obtain the solvability condition 



dTp\ 
drj ) 



{^5V^ 



(TqK 



drj 



drj^ 







which implies 



Ti=0. 



Consequently, ()5.1ip yields 

-ar| (exp 
Integrating over ^, we get 



\,e 



2(71 



973 

5? 



exp 



2ai 



orj 



dT3 _ gio ?7 ^tq 

9e " 



Ao drj ■ 
Using dESJ, ([5T21) and ([5T0| in dO]), we have, at 0(£2/3)^ 



Aoe 



dT4 

~d£: 



fji 






0-2 ?7 
(717 77 



>973 
)9t2 

9r; 
9Vo 



0-3 K 



0-6 K 



dT3 
d^ 
dT2 

drj 



■0-4 



-CTs 



9^3 

d^drj 

d^T2 

drj^ 



2dTo 



drj'^ drj 



cri9 ?/ 



drj ' 



(5.10) 



(5.11) 



(5.12) 



(5.13) 



(5.14) 
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Using (|5.13p on the right hand side to ehminate the term dr^/d^ and 
d'^TQ/drf, we get 



to ehminate 



d 
-ai— [exp 



2cri 



dT4 



= exp 



2cri 



10-5?/ +aeK 



9t2 
drj 



c^s 



9S 



I T ^2 '''18 O'l \ t^''o 



Ao J drj j' 



where 



UJ2 



0-4 O-IO , cri7 

Wo = -: 1 : 

Ao erg (T8 

CT3 CTio 0-4 (Tg fJio crgcriv 



CJl 



(^4 gjO + q"l8 gl 
Ao 



(5.15) 



wa 



0'2 CTlO , cr4 fS f 10 f 5 cri7 

cTig 



Ao AoCTs CTs Ao AocTs 

Integrating over ^ in (— cxd, cxd), we obtain the solvabihty condition 



0-8 



(0-5 ry^ + 0-6 K^ 



Wo ?7 - (Wl + ^2 K 77 + W3 77 j 



- as- 



(5.16) 



977 ' ° 977^ "' ^'^' ^ '' '^ ' ' drj 

This is the equation which determines T2, the first non-zero correction term to tq. Let 
us now solve ()5.9p and (I5.16P for To and T2. Since (15. 9p is of the same form as the 
equation p.36p . we can use the same technique as in Section 13.31 to analyse it. We 
find (compare with (|3.38|) ^ 



lim To(r/) = 7^3^/^ 



r;^ — 00 



crfCTB 



1/3 



H\ 



12 



1/3 



(J(, K, 



(5.17) 



where the function i/ : R ^ (0, 00) is given by p.40p . Moreover, we have the following 
identity (compare with (13.371) 1 



exp 



0-577^ + ScTe K77 



3cr8 



V 



dro 1 , 

-TT- = / exp 

077 0-8 



Scts 



iz. 



(5.18) 



Multiplying equation (|5.16p by exp [(cr577'^ + 3(TgK7;)/(3cr8)] , using (|5.18p . integrating 
over 77 and using integration by parts, we obtain 



^- r (^y + ^iv-y) + ^^iv'-y') + i^iv^-y^) 



drj 



UJ2 



X exp 



2 v/ ., ' ^^2-v, y J + 4 2 



0-5(2/^ - f7^) + 3cr6K(7/ - 77) 



3(78 



dy. 



Integrating over 77 in (—cxd, cxd), we get 



lim T2 (7;) 

ri — > — CO 



LOq 



0-8 J- 



2/ exp 



00 ^ —00 



(7^{y^ - z^) +?,aeK{y~ z) 



+ 



(y - z) exp 



00 J —00 

■00 



wi 

W2 

2^ 

^^^8 -'-00 J -00 



3(78 
cr5(y3 - z'^) +i(jQK{y- z 



(y^ -2:^)exp 



00 ^ — 00 



3(78 
(75 (y^ - 2;^) + 3(76k(7/ - z) 



3(78 
0'5(7/'^ - Z^) + 3(76k(77 - Z) 
3(78 



dydz 
dy dz 

dy dz 
dydz. (5.19) 
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Using the substitution 

we simplify the integrals on the right hand side of (j5.19p using the same technique 
as in (|3.38p . In all cases, we are able to integrate over u variable explicitly because 
the u-integral is the Gaussian integral. The last two integrals are zero because we 
integrate odd functions in u over (—00,00). Thus (j5.19p becomes 



hm T2[ri) = [ujo / V^exp 



1/3 
-V" — I 7T I (Tq KV 



.3 ' 12 



dv. (5.21) 



Let us define the function G : M ^ (0, 00) by 

/•oo 

G{z)^ I Vvexp[-'y3 + zw] dw. (5.22) 

"'0 

Using (|5.15p and the definition (|5.22p . the formula (|5.2ip implies 



lim r,(,)= ...-^-^^p^G -^ a,A. (5.23) 

The limits (|5.17p and (|5.23p give us the desired estimates of the period of oscillation. 
Transforming back to the original variables, the zero-order approximation to the mean 
period of oscillation is given by 

To(fcO = V^ 3Va (_^) '" H [- [-^y a, (k, _ A')) (5.24) 

where the constants a^ , cg and ag are given in Appendix [C) More precisely, the 
formulae for the constants CT5, erg and trg are obtained by dropping the overbars in 
their definitions in Appendix[C] They only depend on the derivatives of the coefficients 
(|4.2p - (|4.6p and on the eigenvectors Ui and U2. The approximation To(fci) is plotted in 
Figure [?nT a) as the black dashed line. Transforming (|5.23p into original variables, we 
can obtained an improved (second-order) approximation of the period of oscillation 
asTo(fci)+T2(fci) where 

Uh)=(a,,~^-^~'-^)^G(-(^Y\uh-K)] (5.25) 
V Ao Ao / crsCTs \ \a5a^J J 

and the constants cr^, i — 1,4,5,6,8,10,17,18, are given in Appendix [Cl Again, 
the overbars have to be dropped. The approximation To(fci) -I- T'2(fci) is plotted in 
Figure lOT a) as the black solid line. 

In Figure I2.4r b) , we studied the period of oscillation as a function of the volume 
V for kid ~ Kd- To approximate this dependence, we put kid = K into formulae 
(15:2411 and ^2^ . Since G(0) = ^/1^/i and iJ(0) = r(l/6)/3, we obtain 

ro(A-) = V^ 3-5/6 f^\ p(i/g)^ (5 26) 

Ti{K)=\gyj -=, (5.27) 

V -Ao ^Q J asa^\/6 
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Fig. 5.1. (a) Approximations of the mean period of oscillation To{ki) and To(fci) + T2(fci) 
given by \5.24\ and \5.25\ . We also replot the results of Figure \2.4}{ a) for comparison. Parameters 
are the same as in Figure W^a). (b) Approximations of the mean period of oscillation Tq{K) 
and To(K) + T2(K) given by i5.26l and \5.2T\ . We also replot the results of Figure \2^ b) for 



comparison. Parameters are given by i2.5\ 



where F is the Gamma function. The values oiTo{K) and To{K) +T2{K) are plotted 
as functions of the volume V in Figure [SHT b). We see that To{K) + T2{K) compare 
well with the results obtained by the stochastic simulation for kid = Kd- To be 
more precise, we should compare it with stochastic results obtained for kid = K 
where K depends on the volume as shown in Figure I4.4r a') . Such data are plotted in 
Figure [5?lT b) too. 

Another way to approximate the period of oscillation is to approximate dx{x,y, ki) 
(resp. dxy{x,y,ki) and dy{x,y,ki)) by d^iq) (resp. d^yiq) and dy{q)) only in the 
equation (|5.ip . If the approximation of Vx {x, y, fci ) and Vy (a;, y, fci ) is the same as 
before, the resulting formulae are equal to (|5.24p ~ (|5.27p with an = 0. The comparison 
of the results for an = is shown in Figure 15.21 We see that we have an excellent 
comparison with the results of the stochastic simulation for kid — K as & function of 
the volume V - panel Figure I5.2r b) . 



6. Conclusion. Bifurcation diagrams of ODEs computed by standard numer- 
ical continuation methods constitute a systematic way to study and summarize the 
dependence of the behaviour of the ODE system on its parameters. An exploration 
of such a dependence could be in principle done by a direct integration in time, but 
it would be much more computationally intensive than the numerical continuation. 
In a similar way, the exploration of the dependence of the behaviour of a stochastic 
chemical model on its parameters can in principle be done by the Gillespie SSA; yet, it 
is more efficient to study it by solving the underlying stationary Fokker-Planck equa- 
tion numerically or by analysing it. For example, computation of the mean period 
of oscillation by the Gillespie SSA is more computationally intensive for large values 
of the system size (volume V) because the large values of V correspond to chemi- 
cal systems with many molecules. On the other hand, the computational intensity 
of the FEM solution of the Fokker-Planck equation is independent of V, making it 
more suitable for studying the dependence of the mean period of oscillation on the 
volume V . The asymptotic formulae (|5.24p and (|5.25p provide further insights into 
the behaviour of the mean period of oscillation. At the bifurcation point kid — K 
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Fig. 5.2. (a) Approximations of the mean period of oscillation To(fei) and To{ki) + T2{ki) 
given by ^5.24i and i5.25t with an = 0. We also replot the results of Figure \2.4^ a) for comparison. 
Parameters are the same as in FigureTK^a). (b) Approximations of the mean period of oscillation 
To{K) and To{K) + T2{K) given by i5.26l and ]5.27l with an = 0. We also plot the results obtained 
by the stochastic simulation for fci^ = K. The other parameters are given by 112.51 1. 



of (|4.17p . the formula (|5.24p simplifies to (|5.26p . Inspecting (|5.26p . we conclude that 
the mean period of oscillation asymptotes to infinity as V^^^ for V — > oo. If V is 
finite, the stochastic model of the chemical system (|2.ip - (l2.2p possesses oscillatory 
solutions for any value of kid which is close to K or close to Kd where kid = Kd is 
the bifurcation point of the classical deterministic ODE model (|2.3p - (l2.4p . On the 
other hand, the trajectories of (|2.3p - (|2.4p only oscillate for kid > Kd- The period of 
oscillation asymptotes to infinity as {kid — Kd)~^^'^ for kid -^ K^ . 

In this paper, we used an illustrative example of the chemical system with two 
chemical species for which the deterministic mean-field description is undergoing a 
SNIPER bifurcation. The system was simple enough that we could directly compute 
many realisations of the Gillespie SSA. Averaging over many realisations, we estimated 
important characteristics of the system, for example, the mean period of oscillation. 
These estimates were used for the visual comparison with the results obtained by the 
asymptotic analysis of the chemical Fokker-Planck equation and by solving this equa- 
tion numerically. Another advantage of the illustrative chemical system was that the 
chemical Fokker-Planck equation was two-dimensional, which simplified its asymp- 
totic analysis and numerical solution. If we have a system of N chemical species, the 
resulting chemical Fokker-Planck equation will be iV-dimensional. We are currently 
investigating the advantages and disadvantages of this approach for N larger than 2. 
Clearly, a suitable numerical method for solving the higher-dimensional PDEs must 
be applied. The analysis presented here can also be generalized to higher-dimensional 
cases. All we need to find are the eigenvectors of the Jacobian matrix at the saddle. 
The biggest contribution to the mean exit time or the mean period of oscillation will 
be in the direction of the eigenvector corresponding to the saddle-node connection. 
Thus the approach presented is potentially equally applicable to the A^-dimensional 
case and lead to the estimates of the dependence of the important system charac- 
teristics (e.g. the mean period of oscillation) on the model parameters. Alternative 
methods to obtain these estimates include: accelerating the Gillespie SSA by us- 
ing approximate SSAs [S], or estimating the low-dimensional effective Fokker-Planck 
equation by using short bursts of appropriately initialized stochastic simulations [3] . 
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Both approaches use SSAs. In contrast, the methods presented in this paper were 
based on the numerical solution and asymptotic analysis of PDEs, i.e. no stochastic 
simulation and no generators of random numbers were needed. 

Appendix A. Chemical Fokker-Planck equation for TV-dimensional sys- 
tems. Let us consider a well-stirred mixture of A^ > 1 molecular species that chem- 
ically interact through M > 1 chemical reactions Rj, j — 1, . . . ,M. The state of 
this system is described by X = [Xi, . . . , Ajy] where Xi = Xi{t) is the number of 
molecules of the i-th chemical species, i = 1, . . . , A. Let q;j(x) be the propensity func- 
tion of the chemical reaction Rj, j = 1, . • . , M, i.e. aj{x)dt is the probability that, 
given X(t) = x, one Rj reaction will occur in the next infinitesimal time interval 
[t,t + dt). Let i/ji be the change in the number of Xi produced by one Rj reaction. 
Then the chemical system can be described by the chemical master equation 



M 



d r 

— p(x, t) = j2 ["j(x - i^j)p{^ - I'j.t) 



aj(x)p(x,i) 



where p(x, t) is the probability that X(i) = x and Uj = [vji , . . . , i^jn] ■ Considering the 
integer valued vector x as a real variable, Gillespie [7j derived the chemical Langevin 
equation by approximating Poisson random variables by normal random variables. 
This approximation is possible if many reaction events happen before the propensity 
functions change significantly their values, see [7] for details. The chemical Langevin 
equation can be written in the form 



M 



M 



dX, = K]^,,a,(X(t)) dt + J2'^,^o^y'int))dW, 

\J = 1 / .7=1 

which corresponds to the chemical Fokker-Planck equation 



N 






dt 



E 



OXiXff; 



M 



J^J^jiajix) P(x,i) 



VJ = 1 



1 d^ 

2dx? 



M 



^v,iV,ka,{x) ) p(x,t) 

\3 = ^ 



M 



^^%OLj{x) ] p(x,t) 
^J = l 



(A.l) 



Appendix B. Equation for the exit time r. Let p{x' ,y',t;x,y,0) be the 
probability that X{t) = x' , Y{t) — y' given that A(0) = x and 1^(0) = y. It satisfies 
the backward Kolmogorov equation 



dp d'^p d'^p d'^p dp dp 

2; T^ 9 i CLxy "7^ 7^ "r O^y "TT ^ -r Vx tt -p Vy - 



dt 



dx^ dxdy dy"^ dx dy 



(B.l) 



Let h{x,y,t) be the probability that [X{t),Y{t)] € 17 at time t given that it started 
at [A(0),r(0)] = [x,y]. Then 

Hx,y,t)^ / p{x',y',t;x,y,0)dx'dy'. 



Integrating (|B.1[) . we obtain 



dh d^h d^h d^h dh dh 

ot ox'^ oxoy ay^ ox oy 



(B.2) 



STOCHASTIC BIFURCATIONS 31 

The mean exit time r to leave fi, given that initially [X(0),F(0)] = [x,y], can be 
computed as follows [5] 

f°° dh f°° 

T{x,y)^~j t—{x,y,t)dt^ h{x,y,t)dt. 

Integrating equation (jB.2p over t, we obtain the following elliptic problem: 
d'^T d'^T d'^T dr dr r r i r^ 

Appendix C. Formulae for coefficients q and Ui. The coefficients Ci{dx): 
Ci{dxy), Ci{dy), CiilJx) and Ci(wj,) can be computed by the Taylor expansion. Below, 
we provide the expressions for those coefficients which actually appear in the final 
formulae for the mean period of oscillation: 






V d V 

csivy) = ^Tzir (?) -;r' + o-^ (g)"i2M22, ci{dx) = ^^(g)uii + ^T^(g)u21, 



d^Vy u?2 , 9^Vy _ - ddx ,_, _ , ddx ,_, _ 



ci[dxy) = ^j^(9)uii + -7^{q)U2i, ci[dy) = -^(gjun + ^p(g)u2i, 

where we used the formulae (|4.2|) ~ (|4.3|) to simplify the Taylor expansion oivx{x, y, ki) 
and Vy (x, y,ki). Note that some of the second derivatives of Vx and Vy are zero which 
makes the resulting formulae shorter. However, even if they were nonzero, the same 
analysis could be carried through. The a coefficients are given by 

CTi = {detuy"^ (dxi^)u2i -dxy{'q)u2iuii +dy(g)un) , 

0-4 = (dct[/)"^ {-dxiq)2u2lU22 +dxy{q) (U11U22 + U12U21) -dy{q)2uilUl2) , 
0-5 = {detUy^ {ci{Vx)u22 - Ci{Vy)ui2) , 

(76 = -~{detuy^ui2, 

as = {detUy"^ (dxiq)ul2 -dxy{q)ui2U22 + '^y (9)^12) I 

O-IO = (det C/)^"^ (C2(wa;)u22 - C2(lJy)ui2) , 

0-17 = (det Uy"^ {ci(dx)iq)ul2 ~ ci(dj;,j)(9)Mi2U22 + ci{dy){q)ul2) , 

0-18 = (det Uy^ {C3{vx)u22 - C3(TJy)ui2) . 

To get the coefficients at in the final formulae ()5.24p - ()5.27|) . we simply drop the 
overbars in the expressions above. 
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